;+
;PURPOSE
;	make plots and pull out histogram data from a SLUG data file 
;	This program works with the outputs histograms NOT the files
;	that contain all the clusters
;SYNTAX
;	slug_histogram, filename[, psname=psname, outname=outname, 
;			thresh=thresh]
;INPUTS
;	filename: name of file output by read_slug.pro
;	psname: name of ps file to be written [default to 'slug_histo.ps']
;	outname: name of output structure with info about the histograms
;	ind: which histogram do you want to make from the structure file
;
;Written by R. da Silva, UCSC, 3-16-10
;-
pro slug_histogram, ind, filename , psname=psname, outname=outname, over=over
if n_elements(ind) EQ 0 then ind = 8
if not keyword_Set(filename) then begin
print, 'no filename set using default... BE CAREFUL!!'
filename='slug_histo.fits'
endif
if not keyword_set(outname) then outname='slug_histo_struct.fits'
a=mrdfits(filename,1)
if not keyword_set(psname) then psname='slug_histo.ps'
!p.multi=[0,2 , 2]
!p.font=0
!p.charsize=2
ntime=n_elements(a[0].time_range)
nquan=n_elements(a[0].quan_mid)
if not keyword_set(over) then nquan-=2
psopen, psname, /landscape

tbl=a[ind].tablebins[1:nquan, *]
wh=where(total(tbl, 2))
min1=min(wh)
max1=max(wh)

for i=1, ntime-1 do begin

y=a[ind].tablebins[1:nquan, i]
x=a[ind].quan_mid[1:nquan]
binsiz=x[2]-x[1]
y=float(y)/total(y)
x=x[min1:max1]
y=y[min1:max1]


plot, x, y, ps=10, xtitle=a[ind].histoname, $
	title=rstring(a[ind].time_range[i]/1d6,format='(f6.2)')+' Myr'
for ii=0, n_elements(x)-1 do begin
    oplot, [x[ii]-binsiz/2., x[ii]-binsiz/2.], [0, y[ii]]
    oplot, [x[ii]+binsiz/2., x[ii]+binsiz/2.], [0, y[ii]]
endfor

endfor
!p.multi=[0, 2, 2]
for iii=0, ntime/10 -1 do begin
timesteps=a[ind].time_range
c1=round(ntime/10)*iii
c2=ntime-1
if not keyword_set(over) then c3=1 else c3=0

bbb=a[ind].tablebins[c3:nquan, c1:c2]
bb=(total(bbb, 1) > 1)
b=bb##replicate(1, 1+nquan-c3)
bbb/=b

mea=total(bbb, 2)/(1+c2-c1)
err=stddev_2(bbb)
mea=mea[min1:max1]
err=err[min1:max1]

ploterror, x, mea, replicate(0, nquan), $
                err, ps=10, $
      title=rstring(timesteps[c1]/1d6, format='(f6.1)')+' to '+$
        rstring(timesteps[c2]/1d6, format='(f6.1)')+' Myr'
endfor

psclose




end
